Dynamics of non-equilibrium thermal entanglement 
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The dynamics of a simple spin chain (2 spins) coupled to bosonic baths at different temperatures 
is studied. The analytical solution for the reduced density matrix of the system is found. The 
dynamics and temperature dependence of spin-spin entanglement is analyzed. It is shown that the 
system converges to a steady-state. If the energy levels of the two spins are different, the steady-state 
concurrence assumes its maximum at unequal bath temperatures. It is found that a difference in 
local energy levels can make the steady-state entanglement more stable against high temperatures. 
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I. INTRODUCTION 

In describing real physical systems one should always 
take into account the influence of the surroundings. The 
study of open systems is particularly important for un- 
derstanding processes in quantum physics [i|. Whereas 
in most cases the interaction with an environment de- 
stroys quantum correlations within the system, it is well 
known that in some situations it can also build up en- 
tanglement I4I and in principle even prepare complex 
entangled states [3|- The dynamics of entanglement in 
open systems provides many interesting insights into re- 
laxation and transport situations, in particular if the sys- 
tem dynamics involves many-body interactions (such as 
spin chains, see j4| for a review). In order to understand 
the role of the various parameters that compete in this 
setup, it is useful to find exactly solvable models. Here 
we study the dynamics of a model that was recently in- 
troduced by L. Quiroga @. It consists of a simple spin 
chain in contact with two reservoirs at different temper- 
atures. In such a non-equilibrium case most studies are 
restricted to the steady-state to which the system con- 
verges in the limit of long times [a, H, 0, Q ■ The dynamics 
for the model in the zero-temperature limit was studied 
in [9|. In the following, we study the dynamics of this 
model for generic temperatures. 

This paper is organized as follows. In Sec. II we de- 
scribe the model of a spin chain coupled to bosonic baths 
at different temperatures as introduced in Ref. [5|. For 
completeness we follow [5] in deriving the master equa- 
tion for the reduced density matrix in the Born-Markov 
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approximation. In Sec. Ill we present the analytical so- 
lution for the system dynamics and show the convergence 
of the obtained solution to the density matrix of the non- 
equilibrium steady state solution. One should note that 
in [^ this steady-state was found only in the case when 
the energy levels of the spins are equal. Finally, in Sec. 
IV we discuss the results and conclude. 



II. MODEL 

We consider the simplest spin chain consisting of two 
spins, with each spin coupled to a separate bosonic bath. 
In the derivation of the master equation we follow the 
formalism suggested in Ref. Q. The total Hamiltonian 
is given by 

H = Hs + Hbi + Hb2 + HsBl + HsB2, 



where 



H, = '-la 



£2 



—al + K{(Tja^ + (Ti al 



is the Hamiltonian describing spin-to-spin interactions 
and al , a^ are the Pauli matrices. In this paper units are 
chosen such that kB — h — 1. The constants ei and 62 
denote the energy of spins 1 and 2, respectively, whereas 
K denotes the strength of spin-spin interaction. We will 
see later that the energy difference Ae — ei — £2 has a cru- 
cial role in determining the entanglement of the thermal 
state. We refer to the case Ae = studied in Q as the 
symmetric case. Our study focuses on the non-symmetric 
case Ae 7^ 0. The Hamiltonians of the reservoirs for each 
spin j = 1,2 are given by 
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Hbj - 



E' 



'"j"ri,j"n,. 



The interaction between the spin subsystem and the 
bosonic baths is described by 



The operators Vj^f^ are chosen to satisfy [Hs,Vj,fj,] = 
ujj^^Vj^fj,, and the /j.^ act on bath degrees of freedom 
(this is always possible; their explicit form will be given 
later on). Physically, the index fj. corresponds to transi- 
tions between eigenstates of the system induced by the 
bath. The whole system (spin chain with reservoirs) is 
described by the Liouville equation 



where k — \ K"^ 



dt 



-i[H,a]. 



We assume that the evolution of the dynamical subsys- 
tem (coupled spins) does not influence the state of the 
environment (bosonic reservoirs) so that the density op- 
erator of the whole system a{t) can be written as: 

a(t)=p(t)Si (0)^2(0) 

(irreversibility hypothesis), where each bosonic bath 
is described by a canonical density matrix Bj — 
g-/3jffBj^^j,jg-,9j//Bjj g^j^j -^-j--^ denotes the reduced den- 
sity matrix of the spin chain. 

In Born-Markov approximation the equation for the 
evolution of the reduced density matrix [10] is: 



dp 
'dt ~ 

with dissipators 



-i[Hs,p]+Ciip)+C2{p) 



fi.y 

-{l-e0^---^)[V,,„Vlp]} 
and where the spectral density is given by 

To find a solution we go to the basis of the eigenvectors 
|Ai) with eigenvalues Xi of the Hamiltonian Hs, 



|Ai) = |0,0), Ai = - 



ei + £2 



|A2) = |1,1), A2 



ei + £2 



|A3)=cos(0/2)|l,O)+sin(0/2)|O,l), A3 = k. 



IA4) = -sin(6'/2)|l, 0) + cos(0/2)|O, 1), A4 = -k, 



i^ and tan6i = 2K/{Ae). In this 
representation the dissipative operator Ci{p) becomes 

2 

^.(p) = E J^'H-^^M'^v^.^pvl - {p.v]J,^,}+)+ 

with transition frequencies 

UJl = A2 — A3, 
<^2 = A2 + A3 

and transition operators 

Vi,i=cos(0/2)(|Ai)(A3| + |A4)(A2|), 

Vi,2=sin(0/2)(|A3)(A2|-|Ai)(A4|), 
l>2,i=sin(0/2)(|Ai>(A3|-|A4>(A2|), 

1^2,2 =COs(^/2)(|A3)(A2| + |Ai)(A4|). 

In this paper we consider the bosonic bath as an infinite 
set of harmonic oscillators, so the spectral density has the 
form J^^\uj^) = 7j(w^)nj(w^), where nj(w^) = (e'^^^f - 
1)"-^ and J^J^—Ljfj^) — e'^j"^' j(-'^(w^). For simplicity we 
choose the coupling constant to be frequency independent 
71 (w) = 71 and 72(w) = 72- In the basis |Ai) the equation 
for the diagonal elements of the reduced density matrix 
is given by 



Pii(t) 

_^ I P22(t) 
dt I P33(t) 
/944(t) 



B 



i Pii{t) 

P22(i) 

P33(0 

\ P44(i) 



where i? is a 4 x 4 matrix with constant coefficients. The 
time-dependence for the non-diagonal elements has the 
following form 



P^jW=e*'*'^>^j(0), 

where Sij is a complex number. For the initial state of 
system in the computational basis {|00), |01), |10), |11)} 
we choose 

p(0) =Po|00)(00| +pi|01)(01| +p2|10)(10| 

+ (1 - po - Pi - P2)|ll)(ll| + Ci2|01)(10| + cj2|10)(01|. 



III. EXACT SOLUTION 

The analytical solution in the basis of eigenvectors |Ai 
is given by: 



Pit{t) 



an ai2 aid ou \ / pii(O) 

J I 0,21 0.22 a23 a24 I I P22(0) 

X1Y2 \ 031 032 033 034 I I P33(0) 

O41 042 043 044 / V P44(0) 



where the coefficients a.y are given by: 

an = {X+ + X^e-'''^){Y+ + K,-e-*^^), 



Pisit) = /O34 = e 



i2t\3 



_t(Xl±Z2l 



P43(0). 



In the present solution we have introduced some con- 
stants: 

X- = X+ + X-, 



Y,^Y^+Yr 



X^ = 2cos2(6l/2)j(i)(±a;,) + 2sm^ (9/2) J '■^\±LOi) 



Y^ = 2sin2(0/2)j(i)(±c^,) + 2cos2(0/2)j(2)(±c^,) 



a^2 = {l-e-'^^)il^e-''^)X^Y,+ , 



-tY2) 



or 



tl2^ 



013 = (1 - e-'^^)X^{Y,+ +Y,-e-"^), 



ai4 = {Xt + Xre-*^i)(l - e-'''')Y+, 



a2i = il-e-'^^)il-e-"^)X^Y 



tY2\ 



0.22 



{X-+Xte-'^^){Y^-+Y+e-'^% 



a23 = (X- + X+e-*^0(l - e'''')Y2-, 



024 = (1 - e-*^O^r(V + Y2+e-*''-), 
031 = (1 - e-*^0^r(>^2+ + >^2"e"*''^), 



a32 = (^r + ^i+e-*''^!! - e-*^^)r2+, 



a33 = (Xf + X+e-*^0(>^2+ + Y2-e-'^^), 
a34 = (l-e-*^0(l-e-*^^)Xrr2+, 



041 = {X+ + Xf e-*^0(l - e" ')i^2", 



042 



il~e-'''^)X+iY2-+Y+e-'^^), 



Xf = (jW(±u;0 + ^^'H±w.)) 



^4i^2 + (^g)2 



(J(i)(±c.,)-J^'^(±t^.)) 



^4ii:2 + (Ae)2 

One can easily see that this solution converges with in- 
creasing time to a diagonal density matrix which does 
not depend on the initial conditions: 



lim pu{t) = — — 

t^co X1Y2 



( XtY,^ \ 

^lY2 

X^Y+ 
V X+Y,- J 



lim p-i4.{t) = 0. 

t — ♦CXD 

In the symmetric case (Ae = 0) the above limit in repro- 
duces the result obtained by Quiroga in [5|. In order to 
quantify the entanglement between the spins we consider 
the concurrence [ll|. In the steady-state {t -^ 00) it is 
given by 

Coo = -^M^^{Q, ^|X+K," - X^Y^+\ 



^1 ^1 Y2 Y2 ) 



a43 = (l-e-*^0(l-e-*^^)Xi+V, 



IV. RESULTS AND DISCUSSION 



Taking into account the initial conditions, the non- 
vanishing non-diagonal elements are: 

P34W = e-*2*^^-^^^p34(0), 



The dynamics of entanglement is analyzed in Figures 
1 and 2. In Figure 1 the dynamics of the concurrence 
between the two qubits is shown. For the model con- 
sidered here the spin chain Hamiltonian Hg can entan- 
gle the qubits for specific times, which gives rise to the 
oscillations of concurrence one observes for short times 
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FIG. 1: Dynamics of the concurrence C{t) for the initial re- 
duced density matrix po — \1, 0)(1, Oj. The parameters of the 
model chosen to be "fi — y2 ~ 0.02, ei = 2, €2 = 1, K = 1 
for different temperatures of baths: curve (1) corresponds to 
Ti = 0,5;T2 = 0,2; curve (2) Ti = 1;T2 = 0,5; curve (3) 
Ti = l,5;r2 = l. 
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FIG. 2; Dynamics of the concurrence C{t) for different ini- 
tial states of the reduced density matrix of qubits; Ti = 1, 
Ta = 0.5, 7 = 0.02, ei = 2, ez = 1, K = 1. The curve (1) 
corresponds to po = (|1,0> - |0, 1))((1,0| - {0, l|)/2; curve 
(2) corresponds to po = |1,0)(1,0|; curve (3) corresponds to 
/5o = |l,l){l,l|. 
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FIG. 3: Steady-state concurrence Cao(TM,AT) as a func- 
tion of the mean bath temperature Tm = (31 + T2)/2 and 
temperature difference AT = Ti — T2 in the symmetric case 
ei = £2 = 3 with K = 1. 




FIG. 4: Steady-state concurrence Coo{Tm,AT) as a function 
of the mean bath temperature Tm ~ (31 + 32)/2 and the 
temperature difference AT — Ti — T2 in the case ei = 3, 
£2 = 1, A" = 1. 
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FIG. 5: Maximally possible steady-state concurrence 



K : 



'(ei,e2) in the strong coupling case ei,e2 < K with 



(note that the initial state is chosen to be separable 
Po = |1,0)(1,0|). For large times, the system converges 
to its steady-state. One can see the disappearance of en- 
tanglement with increasing temperatures of the bosonic 
baths which was shown for the steady-state in [5|. In 
Figure 2 the dynamics of the concurrence for different 
initial states of the qubits in shown. For all cases the 
system converges to one and the same value of entan- 
glement. The plots in Figures 1 and 2 show clearly the 
competition between unitary and dissipative dynamics. 
If the qubits starts from a "symmetric state", i.e. |1, 1), 
no oscillations in the concurrence dynamics are observed 
and if the qubits start from a "non-symmetric state" , i.e. 
|1,0)(1,0| or (|1,0) — 10, l))/v^, one can see oscillations of 
the concurrence which correspond to the energy exchange 
between the qubits in the unitary evolution. Both figures 
reveal that after time of order t '^ 2/7 the concurrence 



CT(^ue2) 




FIG. 6: Maximally possible steady-state concurrence 
CaJ^^i^i, £2) in the weak coupling case ei, £2 > K for K — 1; 
In the corner: profile of the 3D surface at the line £1+62 ~ 4. 



"forgets" about initial conditions and converges to the 
same value, given by Coo from the end of the Sec. III. 
The steady-state concurrence Coo(ei, £2, K, Ti,T2) is an- 
alyzed in Figures 3, 4, 5 and 6. In Figures 3 and 4 we 
plot the steady-state concurrence for the symmetric and 
non-symmetric cases as a function of the mean tempera- 
ture {Tm — (Ti -1- T2)/2) and the temperature difference 
(AT = Ti — r2) of the baths. In the symmetric case one 
can easily see that the maximal value of the entanglement 
is reached for equal bath temperatures (AT — 0) 



C!5„ 



sinh(l/T) - 1 



2cosh(cJi/2r) cosh(w2/2r) 



The critical temperature in units of K above which 
the steady-state becomes separable is given by Tc — 
arcsinh(l)^^ {Tc « 1.136) [^]. It is interesting to note 
that in the non-symmetric case (Figure 4) the maxi- 
mal entanglement is reached in the non-equilibrium case 
(AT 7^ 0). In particular, the maximal entanglement is 
larger than the corresponding non-symmetric equilibrium 
concurrence 



^non — sym 



sin 61 sinh((a;2 - uJi)/2T) - 1 
2cosh(wi/2T)cosh(w2/2T) ' 



The temperature at which entanglement disappears is a 
function of the energy difference Ae between qubits: 



Tf7 = 



VAe2/4 + 1 



It is easy to see that this function reaches its minimum 
value in the symmetric case (Ae = 0). In Figures 5 and 6 
we show the maximally reachable value of entanglement 
as a function of qubits energies in the strong and weak 
coupling case. For every pair of energies (ei, £2) we max- 
imize the value of the concurrence for the different tem- 
peratures of the baths (Ti,T2). One can see that in the 
strong coupling case (ei,e2 < K; Figure 5) the maximal 
value of the entanglement corresponds to the symmetric 
case. In Figure 6 one can see that in the weak coupling 
case (ei, £2 > K) the maximal value of the entanglement 
is reached in the non-symmetric case. 

In conclusion, we have found an analytical solution for 
a simple spin system coupled to bosonic baths at differ- 
ent temperatures. We studied the dynamics of the system 
and showed that on the long term the system converges 
to the steady-state solution. Resolving the entanglement 
dynamics allowed us to distinguish between entanglement 
created by the system and by the bath. For the symmet- 
ric case (ei = £2) we reproduced the steady-state found 
by [5j. We focused on the non-symmetric case (ci j^ 62) 
where we found that the steady-state concurrence as- 
sumes its maximal value for unequal bath temperatures. 
This is corresponds to a dynamical equilibrium, where 
the spin chain transfers heat between the baths. We also 
found that a difference in local energy levels can make 
the steady-state entanglement more stable against high 
temperatures. These analytical results motivate further 
numerical studies on longer spin chains. 



arcsinh ^/Ae^/T+l 
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